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ABSTRACT 

Motivated by the recent discovery of massive planets on wide orbits, we present a mechanism for 
the formation of such planets via disk fragmentation in the embedded phase of star formation. In this 
phase, the forming disk intensively accretes matter from the natal cloud core and undergoes several 
fragmentation episodes. However, most fragments are either destroyed or driven into the innermost 
regions (and probably onto the star) due to angular momentum exchange with spiral arms, leading 
to multiple FU-Ori-like bursts and disk expansion. Fragments that are sufficiently massive and form 
in the late embedded phase (when the disk conditions are less extreme) may open a gap and evolve 
into giant planets on typical orbits of several tens to several hundreds of AU. For this mechanism 
to work, the natal cloud core must have sufficient mass and angular momentum to trigger the burst 
mode and also form extended disks of the order of several hundreds of AU. When mass loading from 
the natal cloud core diminishes and the main fragmentation phase ends, such extended disks undergo 
a transient episode of contraction and density increase, during which they may give birth to a last 
and survivable set of giant planets on wide and relatively stable orbits. 

Subject headings: circumstellar matter — planetary systems — protoplanetary disks — hydrodynamics 
— ISM: clouds — stars: formation 



1. INTRODUCTION 

The likelihood of giant planet formation via direct 
gravitational instability of circumstellar disks around 
solar-type stars has been the subject of intense research 
in the past years. Despite much effort in this field, in- 
creasingly sophisticated numerical hydrodynamics sim- 
ulations and analytical considerations continue to yield 
conflicting results. On one hand, some studies indicate 
that giant planets can form in massive disks, particularly 
in their outer parts where conditions for disk fragmenta- 
tion are less extreme and th e competing core-accretion 
model is less vi able ( e -g.. I Johnson fc Gammid 120031; 
IStamatell os et al.l l2007t iMaver et alTT 2007: Boss 2008; 
[Dodson- Robinson et al. 2009; Nero & Biorkman 2009). 
On the other hand, many studies show that gravitational 
fragmentation is unlikely, particularly in the inner few 
tens of AU due to insufficient disk coolin g; and strong 
stellar /envelop e irradiation (e.g.jMatzner fc Levin 2005; 
iRafikov 2005; "Bolev et al.^ 2006, m vA [Rafikovi 120071 : 
[Stamatellos fc Whitworth 2008; Cai e t al.ll2008f ) 



In spite of a great deal of sophistication, the aforemen- 
tioned studies miss one important aspect — circumstellar 
disks are not isolated in the early embedded phase of 
star formation (hereafter, EPSF). In this stage, they 
are subject to intense mass loading from a natal cloud 
core, which can significantly alter the disk's ability to 
fragment. A self-consistent handling of this process in 
numerical simulations is not easy and requires a much 
larger spatial scale (than just that of the disk). A spa- 
tial resolution of less than 1 AU is usually needed for 
planetary-mass fragments to form and survive. Global 
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numerical hydrodynamics simulations of the gravita- 
tional collapse and fragmentation of molecular clouds 
have demonstrated that forming stars are indeed sur- 
rounded by acc reting ^r avitationally un stable disks (e.g., 
iBate fc Boimei l 2003; Krumholz et al.l 12007). Yet, these 
simulations resolve only massive disks, which, if frag- 
mented, produce brown dwarfs or low-mass stellar com- 
panions rather than giant planets. Moreover, such nu- 
merical simulations are very computationally intensive 
and are unable to explore a wide parameter space and 
long evolution times. 

On the other hand, semi-analytic models and simpli- 
fied numerical simulations of the gravitational collapse 
of dense cloud cores can explore a wide range of initial 
conditions and can give us a valuable insight into the 
required conditions for disk fragmentation. Using the 
thin-disk approximation, we were able to self-consistently 
follow the process of cloud core collapse and star/disk 
formation for at least se veral Myr after the formation 
of a c entral stellar object (|Vorobvov fc Basull2005L I2006L 
l2009l ) . These studies have shown that circumstellar disks 
may be gravitationally unstable and susceptible to frag- 
mentation if the rate of gas deposition onto the disk from 
the cloud core is greater than that from the disk onto 
the star, disk viscosity is not too high, and the natal 
cloud cores are characterized by sufficiently large rota- 
tion rates. More sophisticated numerical hydrodynamics 
simulations, though with an approximate treatment of 
gas infall onto the disk, and semi-analytic studies have 
confirmed the susceptibility of non-isolate d disks to frags;- 
mentation, particularly at large radii (e.g:.,|K ratter et al.| 
2008; Rice at al. 2009; Bolev 2009'; iBolev et al. 2009; 
[Clarkell2009l : lRafikovll2009f ). 

The feasibility of disk fragmentation and giant planet 
formation is only one part of the problem. The other 
part is the likelihood of survival of giant planets formed 
via disk fragmentation. Rapid radial migration due to 
gravitational interaction of a giant planet with a natal 
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gas disk (iGoldreich fc Tremaind fToSOl ) has traditionahy 
been one of the stumbhng blocks for the theory of giant 
planet formation and many mechanisms have been pro- 
posed to stop this mig;ration in the lat e evolution phase 
(see e.g., iTho mmes fc Murrav 2006; Crida fc Mor bidelTil 
120071 : llda fc Li n 2008; Matsu mura et al. 2009). In the 
early EPSF, this problem may be even more severe due 
to the fact that disks are more massive and profoundly 
non-axisymmetric. Indeed, our previous numerical stud- 
ies have shown that fragments forming in the EPSF are 
quickly driven into the inner regions and probably onto 
the protostar due to exchange of angul a r mom entum with 
spiral arms (IVorobvov fc Basul l2005l l2006[ ). We have 
speculated that only those fragments that form in the 
late EPSF, when gravitational instability starts to grad- 
ually decline with time, may have a chance to survive. 

In this paper, we present confirmation that the frag- 
ments formed in the EPSF can survive through this 
extreme phase and form giant protoplanets (hereafter, 
GPPs) on large, relatively stable orbits. This finding 
is made possible by the employment of expanded com- 
putational resources, by improvements in the numerical 
model, and by the use of a wider parameter space in 
comparison to previous works. 

2. MODEL EQUATIONS AND INITIAL CONDITIONS 

Our numerical model is explained in detail in 
IVorobvov fc Basul (|2006f ) and is briefiy summarized be- 
low. We make use of the thin-disk approximation to 
compute the gravitational collapse of rotating, gravita- 
tionally unstable cloud cores. We start our numerical 
integration in the pre-stellar phase, which is character- 
ized by a collapsing starless cloud core with a typical 
radius of 10^ AU, continue into the EPSF, which sees 
the formation of a star /disk/envelope system, and ter- 
minate our simulations in the late T Tauri phase. The 
mass accretion rate onto the disk, Menv, is not a free 
parameter but is self-consistently determined by the dy- 
namics of the gas in the envelope. We introduce a "sink 
cell" at Tgc = 8 AU and impose a free infiow inner bound- 
ary condition. Ninety per cent of the gas that crosses the 
inner boundary is assumed to land onto the central star 
plus the inner axisymmetric disk at r < 8 AU. The other 
10% of the accreted gas is assumed to be carried away 
with protostellar jets. 

The basic equations of mass and momentum transport 
in the thin-disk approximation are 

^ = -V,-{^v,), (1) 

nil 

s^ = -VpP + Sg^ + (v-n),, (2) 

where H is the mass surface density, V = jf^ Pdz is the 
vertically integrated gas pressure , Z is the vertical scale 
height, Vp = VrV -\- V(i)(f) is the velocity in the disk plane, 
9p = Qr^ + 9(i)(p is the gravitational acceleration in the 

disk plane, and Vp = rd / dr -\- <j)r~^ d / d(j) is the gradient 
along the planar coordinates of the disk. The gravita- 
tional acceleration includes the gravity of a central 
point object (when formed), the gravity of the inner in- 
active disk (r < 8 AU), and the self-gravity of a circum- 
stellar disk and envelope. For the kinematic viscosity u 



that enters the viscous stress tensor 11, we use the usual 
a-prescription, with a spatially and temporally uniform 
a = 0.005. The latter ch oice is based on our recent work 
(jVorobvov fc Basul [2OO9I ). i.e., a is small enough to not 
eliminate the burst mode and large enough to drive sig- 
nificant accretion in the late stages of disk evolution. 

Equations ^ and (|2]) are closed with a barotropic 
equation that makes a smooth transi tion from isother- 
mal to adiabatic evolution at S = (|Vorobvov fc Basul 
2006). For the ratio of specific heat s we use 7=1 . 4. Th e 
7=5/3 case was explored in Vorobvov fc Basul (|2006f ). 
The value of Ecr is calculated during the numerical sim- 
ulations as Scr = mn/i ^cr 2Z, where the critical v olume 
number density ncr is set to 10 cm" (Larso^llOOl) 
and the mean molecular weight /i = 2.33. The scale 
height Z is calculated using the assumption of vertical 
hydrostatic equilibrium. We note that Z is an increasing 
function of radius, which makes Scr increase with ra- 
dius as well. In practice, this means that the inner disk 
regions are significantly warmer than the outer regions, 
since the optically thick regime in the inner regions is 
achieved at lower S. This in turn impedes the devel- 
opment of gravitational instability and fragmentation in 
the inner disk, in agreement with more sophisticated nu- 
merical simulations and theoretical predictions that di- 
rectly solve for the energ y bal ance equation (see e.g., 
Stamatellos fc WhitworTB 120081 : [BQlevll2009l : iBolev et al.l 
20091) . The use of a spatially varying is an important 
improvement of the numerical model as compared to our 
previous works. Equations ([1]) and ([2]) are solved in po- 
lar coordinates (r, (j)) on a numerical grid with 256 x 256 
or 512 X 512 (depending on the model) grid zones using 
the method of finite-differences. The radial points are 
logarithmically spaced. 

Initially, cloud cores have surface densities H and an- 
gular velocities Q typical for a collapsing, axisymmetric, 
magnetically supercritical core (|Basulll997t ). with S, cx 
at large radii. Cloud cores are initially isothermal 
and they are characterized by a specific ratio of the rota- 
tional to gravitational energy P = ^rot / 1 ^grav | • We have 
considered the evolution of 82 cloud cores with masses 
Mci = 0.2 - 3.0 M0, energy ratios /3 = 0.2 - 2.2 x 10-^ 
and initial gas temperatures T = 10 — 18 K. The initial 
column density is fiattened near the center and achieves 
an asymptotic large-radius profile T^{r) = kcl/{Gr), 
where k = a/^J/tt. In our previous papers (e.g. IVorobvovl 
I2OO9I ). we took A = 2, but here let it vary in the range 
A=2-8, so that models can be more gravitationally un- 
stable. 

3. FORMATION OF GIANT PLANETS 

3.1. Conditions for fragmentation 

The formation of giant planets via disk fragmentation 
can only take place if the the ratio of the local cooling 
time tc to the local dynam ical time is smaller than 
a few (e.g., G ammid 120011 : iRice at all 120031 : 'Meua et al.l 
I2OO5I ). Our preliminary results with disk cooling, viscous 
and shock heating, and stellar irradiation included in the 
code indicate that this condition is sati sfied in the EPSF, 
at least in the disk's outer regions (| Vorobvov fc Basul 
2010). 

Another requirement for disk fragmentation is the 
Toomre criterion Q = Cs^/{7rGTi) < 1, which states that 
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the gas surface density H should be sufficiently high for 
a disk to fragment. Too high a S, however, may prevent 
fragmentation due to inc reased opacity and cooling time 
(|Nero fc Bjorkmajill2QQ9l ). In other words, there exists 
minimum and maximum values of E between which the 
instability and fragmentation are expected to occur. In 
numerical simulations that form disks self-consistently 
(such as our own), S naturally increases from low to- 
ward higher values during the disk formation phase and 
the disk may pass through the unstable regime. 

The critical density for fragmentation can in principle 
be achieved if the rate of mass deposition onto the disk 
^env is greater than the mass flux in the disk Md so that 
H quickly increases with time (Vorobvov & Basu 2007; 
IVorobvovir2009l : IKratter et al.ll2009l : lBolevl l2009) and this 
phase of intense infall (i.e., the EPSF) lasts for many 
dynamical times (so that the GPPs have enough time 
to form). Simple analytical estimates for a viscous disk 
indicate that M^/Menv ^ 3a/Q for T = 30 K (BolevI 
^009), which implies that the outer disk regions are ex- 
pected to fragment for values of a ~ 10 Numerical 
simulations by IVorobvovl (|2009l ) also show that Menv is 
on average several times greater than Md in the EPSF. 

Observations provide conflicting estimates as to the 
duration of the embedded phase Tem. rang:ing: from 
a few xlO^ yr to a fe w x 10^ yr (|An dre fc Montmerle' 
119941 : lEvans et al.l|2009[ ) . From simple analytical grounds 
it follows that rem should be linearly proportional to the 
initial cloud core mass Md and inverse proportional to 
T^/^ and A. In addition, such effects as rotation and 
magnetic fields may steepen the relationship, particularly 
for cloud cores of solar mass and greater. The numeri- 
cally obtained values for /3 = 1.3 x 10~^, T=10 K, and 
A=2 range from 0.0 3 Myr for Md = 0.2 M© to 0.5 Myr 
for Md = 2.1 Mq (|Vorobvovl 2010). These values may 
decrease by as much as 3-4 for T=20 K or A=8. Never- 
theless, cloud cores in the aforementioned mass range are 
expected to have rem that are considerably longer than 
typical orbital times of 360-3200 yr for radial distances 
of 50-100 AU and stellar masses of 0.1-1 Mq. 

Another important initial parameter that determines 
the disk propensity to fragment is the amount of rota- 
tion in the natal cloud core. Indeed, cloud cores with 
greater rotation rates would make larger disks since the 
centrifugal radius is proportional to the square of the 
specific angular momentum, which would increase disk 
tendency to fragment (recall that fragmentation tends 
to occur at large radii). This has been confirmed by nu- 
merical and semi-analytic simulations of disk formation 
(IVorobvov fc B asu 20Q6|; IKratter et all 120081 : IVorobvovl 
120091 : iRice aral...2009l 

The above analysis indicates that disks formed from 
sufficiently massive cloud cores with sufficiently high an- 
gular momenta are expected to fragment in the embed- 
ded phase, particularly in the outer regions. This con- 
clusion has been corroborated by numerical simulations 
that f orm d i sks self-consiste ntly (Vorobvov fc Basu 2005', 
I2006L 120091 : IVorobvovl 120091 : IKratter et al.t 2009) or im- 
pose some p rescribed rate of mass infall onto the disk 
(jBolevI 120091 : iBolev et al] 120091 : iRice at al.ll2009D . It is 
therefore not so much the feasibility of fragmentation as 
the likelihood of survival of the forming fragments that 
we focus in the present study. 



3.2. Numerical results 

In agreement with Section 13. If we have found that 
the disk propensity to fragment increases along the se- 
quence of increasing Md and /3. This tendency is re- 
flected in both the greater numbers and higher masses 
of the fragments, which form via fragmentation of dense 
spiral arms. The final masses of fragments that form in 
our simulations lie in a wide range starting from a few 
Jupiter masses and ending with low- and intermediate- 
mass brown dwarfs. For disk fragmentation to take 
place, the minimum cloud core mass should be Md, min ^ 
0.8 Mo for /3 = 0.2 X 10"^ and A=2. For /3 = 2 x 10"^ 
and A = 2, the corresponding Md,min is approximately 
0.2 Mq. These minimum masses decrease by about 30% 
for higher perturbation amplitudes A = 8 due to an asso- 
ciated increase in Menv, which promotes disk fragmenta- 
tion. On the other hand, our preliminary results suggest 
that Md,min may increase somewhat if a more accurate 
treatment of the thermal physics is considered. 

Not all fragments evolve ultimately to GPPs, most 
have either migrated through the inner boundary or dis- 
persed in the outer regions (possibly due to insufficient 
resolution). Only 6 out of the 82 models have formed 
GPPs with final masses in the Mpi = 5 — 10 Mj range 
on relatively stable orbits at 25-200 AU. In two models 
we saw two GPPs forming simultaneously, but in both 
cases the outer one had dispersed just after 1.0 Myr of 
evolution, either due to insufficient numerical resolution 
or tidal disruption by the inner GPP because of proxim- 
ity to the 4:1 resonance. These numbers should not be 
treated as representative, since we expect the number of 
survived GPPs to increase in simulations with a higher 
numerical resolution. 

Among 47 models with /3 < 10~^ only two models with 
high density perturbation amplitude A = 8 (as opposed 
to standard A=2) have revealed planet formation. The 
most likely reason why low-/3 and low- A models fail to 
form GPPs is that their disks are too small and the frag- 
ments are forming too close to the star, which lowers 
their chances to survive. In addition, these models have 
lower disk-to-star mass ratios than their high-/3 and high- 
A counterparts, which results in the formation of lower- 
mass fragments. Such fragments need more time to open 
the gap and slow down the fast inward migration in the 
embedded phase. 

Figure [T] presents the gas surface density maps (in 
g cm-2) in model 1 (Md = 0.9 M©, /3 = 1.3 x 10"^ 
A=2) at six evolution times after the formation of the 
disk. The rotation is counterclockwise (note that we 
zoom in at the bottom row). Several fragments condense 
in the outer parts of the spiral arms as early as 0.09 Myr 
after the disk formation, but none of them have survived 
by the end of the embedded phase at ^ 0.16 Myr when 
about 75-80% of the envelope has been accreted by the 
disk. They are all driven into the sink cell via a very 
efficient exchange of angular momentum with the spi- 
ral arms, possibly leading to m ultiple FU Orionis bursts 
(|Vorobvov fc Basull2005Ll2006f ) or forming giant planets 
on very close orbits. The byproduct of these bursts is disk 
expansion due to the conservation of angular momentum. 
When the mass loading from the envelope diminishes and 
the burst phase ends, this expansion is followed by tran- 
sient disk contraction, during which gas surface density 
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Fig. 1. — Gas surface density maps (g cm ^, log units) at six times after the formation of the central star (bright circle in the coordinate 
center) in model 1 (M^i = 0.9 Mq, /S = 1.3 x 10~^, A=2). Note that we zoom in as the time increases. The top two rows contain images 
of size 600 AU on each side, while the bottom row contains images of size 400 AU on each side. 



increases and several more fragments form in the disk's 
outermost regions {t = 0.21 Myr). These fragments 
quickly migrate in the inner 100 AU and, hy t = 0.3 Myr, 
only one fragment survives, which later opens a gap and 
evolves into a well defined GPP possessing its own coun- 
terrotating minidisk. Such counterrotating minidisks are 
seen around many fragments. We believe that this ef- 
fect is caused by the gravitational capture of some of 
the neighboring material, which receives a counterrotat- 
ing twist around the forming fragment due to differential 
rotation of the natal spiral arm. 

Figure [2] shows the GPP's radial position rpi (top), 
mass (middle), and Hill's radius (bottom), as a function 



of time in model 1. Since we do not use sink particles 
for GPPs, these values should be treated as approximate. 
Upon its formation at > 100 AU, the GPP migrates 
quickly in the inner 40 AU and back (this time more 
slowly) to r 80 AU. The outward migration is then 
followed by a gradual inward migration until the GPP 



finally settles at 



52 AU. 



The mass of the GPP is estimated by integrating the 
azimuthally- averaged gas surface density profile around a 
local maximum at the planet location. The middle panel 
of Fig. [2] shows the upper limits on the total mass of the 
inner GPP plus its mini-disk (dashed line) and the mass 
Mpi of the inner GPP (solid line). The latter value is 
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Fig. 2. — Planet's radial position (top), mass (middle), and Hill's 
radius and vertical scale height (bottom), for model 1. 

uncertain to within a factor of unity. The total (planet 
plus disk) mass is always in the giant planet regime, while 
Mpi is around 5±1 Mj. The bottom panel shows the Hill 
radius rn = rpi(Mpi/3(M* + Mpi))!/^ of the GPP (solid 
line) and the azimuthally- averaged vertical scale height 
at the planet position Z (dashed line). It is seen that a 
condition for gap opening, > ^, is satisfied. 

How stable are the orbits of GPPs in our modeling? 
Figure [3] presents the radial position of GPPs vs. time 
in model 2 (solid line, Md = 0.4 M©, /3 = 0.01, A = 2) 
and model 3 (dashed line, Md = 2.0 M©, /3 = 3 x 10"^, 
A = 8). Initially, model 2 demonstrates a similar mi- 
gration pattern to that of model 1. However, the GPP 
in model 2 does not seem to settle at a stable orbit but 
slowly migrates inward. Whether or not this migration 
would stop at later times is unclear. On the contrary, the 
GPP in model 3 appears to stabilize at a rather large ra- 
dial distance of 190 AU. These examples indicate that 
GPPs may have various migration histories, which de- 
pend probably on particular physical conditions in the 



disk and the natal cloud core. The net planet masses are 
5 ± 1 Mj in model 2 and 10 ± 1 Mj in model 3. 
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Fig. 3. — Planet's radial position in model 2 (solid line, = 
0.4 M0, (3 = 0.01, A = 2) and model 3 (dashed line, = 2.0 M©, 
/3 = 3 X 10-3, A = 8). 

We have studied the long-term evolution of disks that 
are formed by the self-consistent collapse of prestellar 
cores. Our model yields gas giant formation starting 
from initial conditions of the early stages of star forma- 
tion. The initial cores are more gravitationally unstable 
and have greater angular momenta than similar models 
studied in the past, and a large number of models have 
been run with relatively high resolution. An early burst 
mode of evolution is characterized by the formation of 
clumps which are then driven into the inner disk. How- 
ever, in a small subset of models, massive fragments are 
formed on wide orbits and settle into stable orbits of ra- 
dius > 50 AU. An interesting feature is that minidisks 
around these fragments can be counterrotating with re- 
spect to the disk. Sometimes, the final orbit can be much 
larger or alternatively there can be a slowly continuing in- 
ward migration. We believe that our results can explain 
the purported observations of giant planets on wide or- 
bits. By extrapolation, they may also represent the first 
stages of the eventual formation of a low mass brown 
dwarf companion. 
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